home *** CD-ROM | disk | FTP | other *** search
/ PC-X 1997 October / pcx14_9710.iso / swag / math.swg / 0093_Simplex Method.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1995-05-26  |  7.5 KB  |  228 lines
  1. {
  2. In response to a request from rfv@maties.sun.ac.za in sunny South
  3. Africa, here is some code from snowy Toronto which implements the
  4. simplex method:
  5.  
  6. {SIMPLEX.PAS -- Simplex routine adapted from FORTRAN version in
  7.                 Numerical Recipes, by Press, Flannery, Teukolsy, and
  8.                 Vetterling, 1986 edition.  Adaptation made by
  9.                 Howard L. Kaplan, Addiction Research Foundation,
  10.                 Toronto, Ontario in 1993.  If you have the right to
  11.                 use the FORTRAN original, then Howard Kaplan grants
  12.                 you the right to use this adaptation.}
  13.  
  14. {$N+}
  15. {if you want a standalone test, then do not $define UnitMode}
  16. {$define UnitMode}
  17.  
  18. {$ifdef UnitMode}
  19. unit    SIMPLEX;
  20.  
  21. interface
  22. {$endif}
  23.  
  24. const   MaxSimplexDimensions=10;
  25.         MaxIterations=200;
  26.  
  27. type    ParameterVector=array[1..MaxSimplexDimensions] of single;
  28.         ParameterArray=array[1..MaxSimplexDimensions+1] of
  29.                        ParameterVector;
  30.         EvaluationArray=array[1..MaxSimplexDimensions+1] of single;
  31.         String20=string[20];
  32.         Evaluator=function(R: ParameterVector;
  33.                            const Why: OpenString): single;
  34.  
  35. procedure Amoeba(var P: ParameterArray; {starting vertices}
  36.                  NDim: integer; {N of elements of ParameterVector used}
  37.                  FTol: single;  {fractional convergence tolerance}
  38.                  Func: Evaluator; {fit evaluation function}
  39.                  var Iter: integer {number of iterations used});
  40. {$ifdef UnitMode}
  41. implementation
  42. {$else}
  43. forward;
  44. {$endif}
  45.  
  46. procedure Amoeba(var P: ParameterArray; {starting vertices}
  47.                  NDim: integer; {N of elements of ParameterVector used}
  48.                  FTol: single;  {fractional convergence tolerance}
  49.                  Func: Evaluator; {fit evaluation function}
  50.                  var Iter: integer {number of iterations used});
  51. const Alpha=1.0;
  52.       Beta=0.5;
  53.       Gamma=2.0;
  54. var   Y:     EvaluationArray;
  55.       PBar:  ParameterVector; {centroid without worst point}
  56.       PR,
  57.       PRR:   ParameterVector; {test points to be evaluated}
  58.       MPts,              {number of points in the simplex}
  59.       I,J,
  60.       IHi,               {index of highest (worst) evaluation}
  61.       INHi,              {index of next-highest evaluation}
  62.       ILo:   integer;    {index of lowest (best) evaluation}
  63.       YPr,               {function evaluated at PR}
  64.       YPrr,              {function evaluated at PRR}
  65.       STemp: single;
  66.       Why:   String[20];
  67. {sub}function NString(N: integer): string20;
  68. var   S: String[20];
  69. begin
  70.      Str(N,S);
  71.      while (Length(S)<4) do
  72.         S:=' '+S;
  73.      NString:=S
  74. end;
  75. begin
  76.      MPts:=NDim+1; {number of points in the simplex}
  77.      Iter:=0;
  78.      FillChar(Y,SizeOf(Y),0); {facilitate debugging}
  79.      for J:=1 to MPts do
  80.         Y[J]:=Func(P[J],'Initial row '+NString(J)); {initial simplex}
  81.      repeat
  82.        {Find the worst, next worst, and best vertices so far}
  83.         if (Y[1]>Y[2]) then begin
  84.            IHi:=1;
  85.            INHi:=2
  86.            end
  87.         else begin
  88.            IHi:=2;
  89.            INHi:=1
  90.            end;
  91.         ILo:=1;
  92.         for I:=1 to MPts do begin
  93.            if (Y[I]<Y[ILo]) then
  94.               ILo:=I;
  95.            if (Y[I]>Y[IHi]) then begin
  96.               INHi:=IHi;
  97.               IHi:=I
  98.               end
  99.            else if (Y[I]>Y[INHi]) then
  100.             if (I<>IHi) then
  101.               INHi:=I
  102.            end;
  103.        {If the worst is 0 or close, return}
  104.         if (Y[IHi]<=FTol) then
  105.            Exit;
  106.        {Compute the fractional range from worst to best, and return if
  107.         satisfactory}
  108.         if ((2*Abs(Y[IHi]-Y[ILo])/(Abs(Y[IHi])+Abs(Y[ILo])))<FTol) then
  109.            Exit;
  110.        {If we are not allowed to do any more work, return although
  111.         unsatisfactory}
  112.         if (Iter=MaxIterations) then
  113.            Exit;
  114.        {Do another iteration}
  115.         Inc(Iter);
  116.        {Compute the centroid of the face that leaves out the
  117.         one worst point}
  118.         for J:=1 to NDim do begin
  119.            STemp:=0;
  120.            for I:=1 to MPts do
  121.             if (I<>IHi) then
  122.               STemp:=STemp+P[I,J];
  123.            PBar[J]:=STemp/NDim
  124.            end;
  125.        {Reflect the simplex from the worst point, and evaluate}
  126.         for J:=1 to NDim do
  127.            PR[J]:=(1+Alpha)*PBar[J]-Alpha*P[IHi,J];
  128.         YPr:=Func(PR,'Reflect worst');
  129.         if (YPr<=Y[ILo]) then begin
  130.           {This is better than the best so far, so try an additional
  131.            extrapolation factor of Gamma}
  132.            for J:=1 to NDim do
  133.               PRR[J]:=Gamma*Pr[J]+(1-Gamma)*PBar[J];
  134.            YPrr:=Func(PRR,'Extend reflection');
  135.            if (YPrr<Y[ILo]) then begin
  136.              {replace the highest point with PRR}
  137.               P[IHi]:=PRR;
  138.               Y[IHi]:=YPrr
  139.               end
  140.            else begin
  141.              {replace the highest point with PR}
  142.               P[IHI]:=PR;
  143.               Y[IHi]:=YPr
  144.               end
  145.            end {YPr<YLo}
  146.         else if (YPr>=Y[INHi]) then begin
  147.           {The new point is worse than the second highest, but it might
  148.            still be better than the highest}
  149.            if (YPr<Y[IHi]) then begin
  150.               P[IHi]:=PR;
  151.               Y[IHi]:=YPr
  152.               end;
  153.           {Whether or not the new point replaced the highest, see
  154.            whether a point of the interior, interpolating between
  155.            the (possibly new) highest point and the old centroid,
  156.            is better than the highest point, and if so, replace the
  157.            highest point; if not, contract around the best point so
  158.            far.}
  159.            for J:=1 to NDim do
  160.               PRR[J]:=Beta*P[IHi,J]+(1-Beta)*PBar[J];
  161.            YPrr:=Func(PRR,'Interp. reflection');
  162.            if (YPrr<Y[IHi]) then begin {replace}
  163.               P[IHi]:=PRR;
  164.               Y[IHi]:=YPrr
  165.               end
  166.            else {contract}
  167.               for I:=1 to MPts do
  168.                if (I<>ILo) then begin
  169.                  for J:=1 to NDim do
  170.                     P[I,J]:=(P[I,J]+P[ILo,J])/2;
  171.                  Y[I]:=Func(P[I],'Contract')
  172.                  end
  173.            end {PRR was worse than the second-highest}
  174.         else begin
  175.           {PR was better than the second-highest, so we replace the old
  176.            highest point}
  177.            P[IHi]:=PR;
  178.            Y[IHi]:=YPr
  179.            end
  180.         until (false)
  181. end;
  182.  
  183. {$ifndef UnitMode}
  184.  
  185. function  TestEvaluation(R: ParameterVector;
  186.                          const S: OpenString): single;
  187. {Try to fit Y=A*Sqr(X)+B*X+C, where A=20, B=3, C=16, for -30<=X<=30}
  188. var   N:     integer;
  189.       Sigma: single;
  190. begin
  191.       Write('Evaluating A=',R[1]:6:3,', B=',R[2]:6:3,', C=',R[3]:6:3);
  192.       Sigma:=0;
  193.       for N:=-30 to 30 do
  194.          Sigma:=Sigma+Sqr((R[1]*N*N+R[2]*N+R[3])-(20*N*N+3*N+16));
  195.       WriteLn(' : ',Sigma:10:2);
  196.       TestEvaluation:=Sigma
  197. end;
  198.  
  199. procedure TestAmoeba;
  200. const Dimensions=3;
  201. var   Simplex: ParameterArray;
  202.       Iter:    integer;
  203. begin
  204.      WriteLn;
  205.      Simplex[1,1]:=1;
  206.      Simplex[1,2]:=1;
  207.      Simplex[1,3]:=1;
  208.      Simplex[2,1]:=-1;
  209.      Simplex[2,2]:=-1;
  210.      Simplex[2,3]:=-1;
  211.      Simplex[3,1]:=1;
  212.      Simplex[3,2]:=3;
  213.      Simplex[3,3]:=5;
  214.      Simplex[4,1]:=6;
  215.      Simplex[4,2]:=4;
  216.      Simplex[4,3]:=7;
  217.      Amoeba(Simplex,3,0.00001,TestEvaluation,Iter);
  218.      WriteLn('Converged to A=',Simplex[1,1]:6:3,', B=',Simplex[1,2]:6:3,
  219.              ', C=',Simplex[1,3]:6:3,' after ',Iter,' iterations')
  220. end;
  221. {$endif}
  222.  
  223. begin
  224. {$ifndef UnitMode}
  225.      TestAmoeba
  226. {$endif}
  227. end.
  228.